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ABSTRACT 

We investigate several key questions of plasma heating in open-field regions of the corona that connect 
to the solar wind. We present results for a model of Alfven-wave-driven turbulence for three typical open 
magnetic field structures: a polar coronal hole, an open flux tube neighboring an equatorial streamer, and an 
open flux tube near a strong-held active region. We compare time-steady, one-dimensional turbulent heating 
models (Cranmer et al. 2007) against fully time-dependent three-dimensional reduced-magnetohydrodynamics 
modeling of BRAID (van Ballegooijen et al. 2011). We hnd that the time-steady results agree well with 
time-averaged results from BRAID. The time-dependence allows us to investigate the variability of the 
magnetic huctuations and of the heating in the corona. The high-frequency tail of the power spectrum of 
huctuations forms a power law whose exponent varies with height, and we discuss the possible physical 
explanation for this behavior. The variability in the heating rate is bursty and nanoflare-like in nature, and we 
analyze the amount of energy lost via dissipative heating in transient events throughout the simulation. The 
average energy in these events is 10^''^' erg, within the “picoflare” range, and many events reach classical 
“nanoflare” energies. We also estimated the multithermal distribution of temperatures that would result from 
the heating-rate variability, and found good agreement with observed widths of coronal differential emission 
measure (DEM) distributions. The results of the modeling presented in this paper provide compelling evidence 
that turbulent heating in the solar atmosphere by Alfven waves accelerates the solar wind in open flux tubes. 


1. INTRODUCTION 

The solar atmosphere is characterized by a temperature pro¬ 
file that has puzzled scientists for decades. The photosphere, 
the visible surface of the Sun, is roughly 6000 K, while the 
diffuse corona above it is heated to millions of degrees. Ob¬ 
servations have not been able to determine the mechanism(s) 
responsible for the extreme coronal heating, though many 
physical processes have been suggested. The same processes 
that might explain the temperature of the corona can also ac¬ 
count for the acceleration of the solar wind. Only a small frac¬ 
tion of the mechanical energy in the Sun’s sub-photospheric 
convection zone needs to be converted to heat in order to 
power the corona. However, it has proved exceedingly diffi¬ 
cult to distinguish between competing theoretical models us¬ 
ing existing observations. Recent summaries of these prob¬ 
lems and controversies have been presented by, e.g., Parnell 
& De Moortel (2012), McIntosh (2012), Klimchuk (2015), 
and Cranmer et al. (2015). 

If solar wind flux tubes are open to interplanetary space, 
and if they remain open on timescales comparable to the 
time it takes plasma to accelerate into the corona, then the 
main sources of energy must be injected at the footpoints 
of the flux tubes. Thus, in wave/turbulence-driven models, 
the convection-driven jostling of the flux-tube footpoints is 
assumed to generate wave-like fluctuations that propagate 
up into the extended corona. These waves (usually Alfven 
waves) are often proposed to partially reflect back down to¬ 
ward the Sun, develop into strong magnetohydrodynamic 
(MHD) turbulence, and dissipate gradually. In this project, 
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we use two models in this regime, ZEPHYR (Cranmer et al. 
2007) and BRAID (van Ballegooijen et al. 2011). These and 
other such models have been shown to naturally produce re¬ 
alistic fast and slow winds with wave amplitudes of the same 
order of magnitude as those observed in the corona and helio¬ 
sphere (see also Hollweg 1986; Matthaeus et al. 1999; Suzuki 
&. Inutsuka 2006; Ofman 2010; Verdini et al. 2010; Chandran 
et al. 2011; Lionello et al. 2014; Matsumoto & Suzuki 2014; 
Woolsey & Cranmer 2014). 

However, between the above ideas and others that claim 
reconnection plays the primary role in heating, nothing has 
been ruled out because (1) we have not yet made the observa¬ 
tions needed to distinguish between the two paradigms, and 
(2) most models still employ free parameters that can be ad¬ 
justed to improve the agreement with existing observational 
constraints. A complete solution must account for all impor¬ 
tant sources of mass and energy into the three-dimensional 
and time-varying corona. Regardless of whether the dominant 
coronal fluctuations are wave-like or reconnection-driven, 
they are generated at small spatial scales in the lower atmo¬ 
sphere and are magnified and “stretched” as they propagate 
up in height. Their impact on the solar wind’s energy budget 
depends crucially on the multi-scale topological structure of 
the Sun’s magnetic held. 

The current generation of time-steady ID turbulence-driven 
solar wind models (e.g. Cranmer et al. 2007; Verdini et al. 
2010; Chandran et al. 2011; Cranmer et al. 2013a; Lionello 
et al. 2014) contain detailed descriptions of many physical 
processes relevant to coronal heating and solar wind acceler¬ 
ation. However, none of them self-consistently simulate the 
actual process of MHD turbulent cascade as a consequence of 
the partial reflections and nonlinear interactions of Alfvenic 
wave packets. Thus, it is important to improve our under¬ 
standing of the flux-tube geometry of the open-field corona 
and how various types of fluctuations interact with those 
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structures by using higher-dimensional and time-dependent 
models to understand the underlying physical processes. 

In Section 2, we describe the modeling used to improve our 
understanding of the geometry of flux tubes in a time-steady 
sense and extending such understanding to time-dependent 
three-dimensional modeling. We focus on three models 
that are typical of common structures in the solar corona 
and present results of each in Section 3 and compare the 
one-dimensional, time-steady results from earlier work using 
ZEPHYR with the new three-dimensional, time-dependent re¬ 
sults. We then analyze the polar coronal hole BRAID model 
further in Section 4. We include a discussion and conclusions 
in Section 5. 

2. MHD TURBULENCE SIMULATION METHODS 
2.1. Previous time-steady modeling 

The magnetic field profiles used for this project are listed 
in Table 1 and are representative of a polar coronal hole, an 
equatorial streamer, and a flux tube neighboring an active re¬ 
gion. The radial profiles of magnetic field strength B{r) for 
these models are shown in Figure la (see also Banaszkiewicz 
et al. 1998; Cranmer et al. 2007). The Alfven travel time from 
the base of the model to a height of 2 solar radii is calcu¬ 
lated using Equation (3), and the expansion factor is defined 
by Wang & Sheeley (1990). The expansion factor is a ratio 
of the magnetic field strength at a height of 1.5 solar radii to 
the field strength at the photospheric base (here, at a height 
of 0.04 solar radii, the size of a supergranule), normalized 
such that an expansion factor of 1 is simple radial expansion, 
and an expansion factor larger than 1 is considered “superra- 
dial” expansion. Key results from these three models using 
ZEPHYR, a one-dimensional time-steady code introduced by 
Cranmer et al. (2007), are presented in Figure 1. 

The polar coronal hole model has the smallest expansion 
from the photosphere to the low corona, and the highest wind 
speed at 1 AU. This relationship was first seen in empirical 
fits to observations (see, e.g., Wang & Sheeley 1990; Arge 
&. Pizzo 2000) and has also been seen in models (see, e.g., 
Woolsey & Cranmer 2014). The equatorial streamer model 
represents an open flux tube directly neighboring the helmet 
streamers seen around the equator at solar minimum, and pro¬ 
duces slower wind at 1 AU. The flux tube neighboring an ac¬ 
tive region has a stronger magnetic field above the transition 
region and an even slower wind speed at 1 AU than the equa¬ 
torial streamer model. In these models, the transition region 
for PCH and EQS is at a height of 0.01 solar radii, and the 
NAR model has a transition region slightly higher, at a height 
of 0.015 solar radii. The density profiles were determined 
self-consistently with the wind speed using mass flux conser¬ 
vation. 

ZEPHYR solves for a steady state solution to the solar wind 
properties generated by a one-dimensional open flux tube. By 
solving the equations of mass, momentum, and energy con¬ 
servation and iterating to a stable solution, the code produces 
solar wind with mean properties that match observations and 
in situ measurements (Cranmer et al. 2007; Woolsey & Cran¬ 
mer 2014). The code’s expression for the turbulent heating is 
a phenomenological cascade rate whose form has been guided 
and validated by several generations of numerical simulations 
and other models of imbalanced, reflection-driven turbulence 
(see, e.g., Hossain et al. 1995; Lithwick et al. 2007; Chandran 
et al. 2009; Oughton et al. 2015). For additional details, see 
Section 3.2 below. 


ZEPHYR can only take us so far, however. It is only by 
modeling the fully 3D spatial and time dependence of the cas¬ 
cade process (together with the intermittent development of 
magnetic islands and current sheets on small scales) that we 
can better understand the way in which the plasma is heated 
by the dissipation of turbulence. We therefore make use of the 
time-dependent modeling of coronal turbulence introduced by 
van Ballegooijen et al. (2011) using the reduced magnetohy¬ 
drodynamics (RMHD) code called BRAID. 

2.2. Including time-dependence and higher dimensions 

Previous numerical simulations of reflection-driven RMHD 
turbulence with the full nonlinear terms include Dmitruk & 
Matthaeus (2003), Perez & Chandran (2013), and prior stud¬ 
ies using BRAID on closed loops (van Ballegooijen et al. 
2011; Asgari-Targhi & van Ballegooijen 2012; Asgari-Targhi 
et al. 2013, 2014). Our version of BRAID uses the three- 
dimensional equations of RMHD (Strauss 1976; Montgomery 
1982; Zank & Matthaeus 1992; Bhattacharjee et al. 1998) to 
solve for the nonlinear reactions between Alfven waves gen¬ 
erated at the single footpoint of an open flux tube. RMHD re¬ 
lies on the assumption that the incompressible magnetic fluc¬ 
tuations (5B in the system are small compared to an overall 
background field Bq. At scales within the turbulence inertial 
range, observations show that (5B is much smaller in ampli¬ 
tude than the strength of the surrounding magnetic field B and 
(5B is perpendicular to B (Matthaeus et al. 1990; Tu & Marsch 
1995; Chen et al. 2012). 

Some implementations of RMHD combine together the im¬ 
plicit assumptions of incompressible fluctuations, high mag¬ 
netic pressure (i.e., plasma /3 <C 1), and a uniform background 
field Bq. However, the flux tubes we model in the upper 
chromosphere and low corona have some regions with /3 « 1 
and a vertical field Bq that declines rapidly with height. It 
is not necessarily the case that RMHD applies in this situa¬ 
tion. However, it was shown in Section 3.1 of van Ballegooi¬ 
jen et al. (2011) that there is a self-consistent small-parameter 
expansion that gives rise to a set of RMHD equations ap¬ 
propriate for the chromosphere and corona. In these equa¬ 
tions, the dominant, first-order fluctuations are transverse and 
incompressible—even when (3 « 1—and gravitational strati¬ 
fication is included to account for the height variation of Bq. 

BRAID uses two cross-sectional dimensions of the flux 
tube and a third dimension along the length of the flux tube, 
which aligns with the background field Bq. Alfven waves are 
generated at the lower boundary by random footpoint motions 
with an rms velocity of 1.5 km s“' and correlation time of 60 
s. These fluctuations are generated by taking a randomized 
white-noise time stream and passing it through a low-pass 
(Gaussian) frequency filter that removes fluctuations shorter 
than the specified correlation time. The time stream is nor¬ 
malized to the desired rms velocity amplitude and split up be¬ 
tween two orthogonal low-kj^ Bessel-function modes of the 
cylindrically symmetric system. The driver modes are shown 
in Figure 2 and are discussed further in Appendix B of van 
Ballegooijen et al. (201 1). 

The magnetic and velocity fluctuations can be approxi¬ 
mated by 

5B = Vj_/!xB, (la) 

,5v = Vi/xB, (lb) 

where B is the full magnetic field vector, whose magnitude 
varies with height, and B is the unit vector along the mag- 
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Table 1 

Models in BRAID 


Modeled structure 

Travel time to 2 Rq 

Expansion factor 

Speed at 1 AU 

Identifier 

Polar Coronal Hole 

770 s 

4.5 

720 km s“* 

PCH 

Equatorial Streamer 

3550 s 

9.1 

480 km s“' 

EQS 

Near Active Region 

4400 s 

41 

450 km s“' 

NAR 



Height obove Sun (Rq) 



Height obove Sun (Rq) 



10 "^ 10 "'‘ 10 "^ 10 "^ 10 "' 10 ° 10 ' 10 ^ 
Height obove Sun (Rq) 



10"^ lO""* 10"^ 10"^ 10"' 10° 10' 10^ 

Height obove Sun (Rq) 


Figure 1. The primary input to ZEPHYR are (a) the magnetic profiles of three representative coronal structures. Results for the bulk solar wind properties from 
each flux tube include: (b) outflow speed, (c) one-fluid temperature, and (d) Alfvenic heating rate. 


netic field. Also, h(r, t) is a height- and time-dependent func¬ 
tion that is analogous to the standard RMHD magnetic flux 
function, and fir, t) is a velocity stream function (sometimes 
called ip in other derivations of RMHD). We can also define 
the magnetic torsion parameter a = -V\h and the parallel 


component of vorticity w = -V^/ (see, e.g., Montgomery 
1982). The functions hir,t) and /(r,f) satisfy the coupled 
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Figure 2. We show here vorticity u} from (a) cosine and (b) sine driver modes 
for the PCH. 

equations 


— +B • (V^w X Vx/) - 
Va [B ■ Va + B ■ (Vxa X Vx/t)]+£>v, (2a) 

|^=B-V/+^+B-(Vx/xVx/t)+A,„ (2b) 

where Hb is the magnetic scale length and D^, and 
describe the effects of viscosity and resistivity (see detailed 
derivation in van Ballegooijen et al. 2011). Both Dy and 
have a hyperdiffusive dependence so that the smallest ed¬ 
dies are damped preferentially and the cascade is allowed to 
proceed without significant damping over most of the k± in¬ 
ertial range. 

In this paper, we have extended previous work using 
BRAID by using an open upper boundary condition instead 
of a closed coronal loop. The model extends to a height of 
Ztop = 2^0. This height was chosen to model as much of the 
solar wind acceleration region as possible, without extending 
into regions where the wind speed becomes an appreciable 
fraction of the Alfven speed, since the RMHD equations of 
BRAID don’t include the outflow speed. To implement the 
open upper boundary condition, we set uj+ (which corresponds 
to the downward Elsasser variable Z+) to zero, whereas in the 
previous coronal loop models that used BRAID, it was set to 
the opposite footpoint’s boundary time stream as described 


above (Asgari-Targhi et al. 2014, and references therein). 

3. TIME-AVERAGED RESULTS FROM THREE MODELED FLUX 
TUBES 

In this paper, we present three models that are representa¬ 
tive of common coronal structures (the same from Table 1). 
Figure 3 gives some of the time-steady background variables 
for the three models, based on the input magnetic field profiles 
shown in Figure 1 . 



Figure 3. Radial profiles for (a) Alfven speed, (b) Alfven travel time, and (c) 
tube radius. 


The Alfven speed, Va = B/^/Anp shows a rapid rise in 
the upper chromosphere, followed by varied behavior in the 
corona depending on the model. The Alfven travel time is 
defined as a monotonically increasing function of height. 


ta(z) = 


r dz' 

Jo Va(z') 


(3) 


where z = 0 is the photospheric lower boundary of each model. 
The BRAID code uses ta as the primary height coordinate. 
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Figure 3 shows the modeled transverse radius of the flux tube, 
which is normalized to 100 km at z = 0 (i.e., a typical length 
scale for an intergranular bright point) and is assumed to re¬ 
main proportional to in accordance with magnetic flux 
conservation. 

Figure 4 provides the time-averaged results from BRAID 
for the PCH. The transition region is at a travel time of 
roughly 210 s and is shown with a dotted line. In Figure 4a, 
we show the magnetic, kinetic, and total energy densities of 
the RMHD fluctuations. We plot these quantities separately 
since no equipartition is assumed. Figure 4b shows the in¬ 
crease in the rms transverse velocity amplitude, Avrms with 
increasing height, which roughly follows the expected sub- 
Alfvenic WKB relation Av^ms cx below the transition re¬ 
gion. The heating rate is also broken up into magnetic (from 
the Dm term), kinetic (from the term), and total for Fig¬ 
ure 4c. Finally, we show the magnetic field fluctuation as a 
function of travel time in Figure 4d. Like the rms velocity, the 
sub-Alfvenic WKB relation cx is followed below 

the transition region but diverges above it. 

The quantities plotted in Figures 4b, 4d, and 5b have gone 
through two “levels” of root mean square averaging. First, we 
take the variance over all kx modes at a given height and time, 
then take the square root. Then, at each height, we take an 
average over the simulation-time dimension using the squares 
of the first level of rms values. It is those quantities that we 
then take the square root of and plot. 

In Figure 5a, we show the magnitude of the Elsasser vari¬ 
ables, Z± = S\ ± SR/^/Airpo, for the PCH model, and show 
the fluctuations Avrms and A/Jj^s in Figure 5b. Note that here 
Abrma is in velocity units, as the actual fluctuations are divided 
by y/Airpo (where po is the background density) for compar¬ 
ison with the velocity fluctuations. Our boundary condition 
enforces that the incoming waves (Z+) have exactly zero am¬ 
plitude at the upper boundary of our model. Figure 5a also 
shows the radial dependence of the Elsasser variables from 
the ZEPHYR model for this coronal hole presented by Cran- 
mer et al. (2007). The ZEPHYR code computes the Alfvenic 
wave energy using a damped wave action conservation equa¬ 
tion that contains the assumption that Z_ ^ Z+. Thus, when 
reporting the magnitudes of the Elsasser variables here for 
direct comparison with the BRAID results, we make use of 
Equation (56) of Cranmer et al. (2007) and correct for non- 
WKB effects by multiplying these quantities by a factor of 
(l + TZ^)/(l-TZ^), where TZ is the reflection coefficient. With 
this correction, there is good agreement between the modeling 
of the PCH for both codes. We use the same correcting factor 
when plotting the Avrms in Figure 5b. Because ZEPHYR as¬ 
sumes equipartition, Avr^s and are equivalent for that 

model. 

Eigure 6 shows the Elsasser variables and the amplitudes 
of magnetic and velocity fluctuations for both the equatorial 
streamer (EQS) and active region (NAR) models to a height 
of z = Q.5Rq, where turbulence has had time to develop in the 
simulation. It is worthy of note that the NAR model shows 
an increase in Z+ above the transition region where the PCH 
and EQS models show a decrease. For these models in Figure 
6b and the results for the PCH (Figure 5b), the shape of the 
Abrma is expected to have a sharp increase at the transition 
region, while the Avrms doesn’t show it (see, e.g.. Figure 9 of 
Cranmer & van Ballegooijen 2005). 

3.1. Energy partitioning in the corona 


In our previous modeling using ZEPHYR (see Section 2.1), 
we assumed equipartition between kinetic and magnetic po¬ 
tential energy densities. With BRAID, we are able to inves¬ 
tigate how far the PCH model differs from this simplifying 
assumption. Eigure 7 shows the ratio of the time-averaged 
magnetic and kinetic energy densities. Above the transition 
region, equipartition is a valid assumption. However, at the 
base of the flux tube, magnetic potential energy dominates, 
and this transitions to a stronger dominance of the kinetic en¬ 
ergy right up to the transition region. 

Also plotted in Figure 7 are six curves showing predic¬ 
tions from non-WKB reflection with a range of frequencies 
between 0.7 and 4.0 mHz. While the general shape is con¬ 
sistent with the BRAID results, it is interesting to note that 
the linear method of computing non-WKB reflection does 
have its limits. These predictions were made using the coro¬ 
nal hole model from Cranmer & van Ballegooijen (2005) as 
a basis. However, that model had a lower transition region 
(ztr ~ O.OO 3 R 0 ), so we have multiplied the height coordinate 
by a factor of 3.4 to match with the ZEPHYR/BRAID coronal 
hole model used in this project. 

Non-WKB theory predicts the kinetic energy to dominate 
in the chromosphere (e.g., Em/E^ < 1 at low heights). This 
was discussed in Appendix A of Cranmer & van Ballegooi¬ 
jen (2005) as an after-effect of the transition from kink-mode 
MHD waves in the photosphere to volume-filling Alfven 
waves in the upper chromosphere. At low frequencies, the 
kink-mode waves are partially evanescent, with two possible 
solutions for the height-increase of Avrms- One of the solu¬ 
tions has Em/Ek < 1 and the other has Em/Ek > 1. However, 
only the solution with Em/Ek < 1 has a physically realistic 
energy density profile (exponentially decaying with increas¬ 
ing height), and this solution also corresponds to a net upward 
phase speed (see also Wang et al. 1995). 


3.2. bleating rates: comparison with time-steady modeling 


We compare the time-averaged heating rates from BRAID 
with the results from the time-steady modeling using 
ZEPHYR. Figure 8a shows the radial dependence of the heat¬ 
ing rate Q, and Figure 8b shows the ratio between the nu¬ 
merically computed heating rates Q with the phenomenologi¬ 
cal heating rate Qphen, as well as comparisons between Q and 
phenomenological heating rates with added correction factors 
(2z.07 and Qioop), described in the following paragraphs. The 
heating rate Qphen is based on the result of a long series of 
turbulence simulations and models (Dobrowolny et al. 1980; 
Grappin et al. 1983; Hossain et al. 1995; Matthaeus et al. 
1999; Dmitruk et al. 2001, 2002; Chandran et al. 2009). The 
analytical expression for this base phenomenological rate is 
given by: 


Qphen — Po 


4Lx 


(4) 


where po is the solar wind density, Z_ and Z+ are the El¬ 
sasser variable amplitudes that represent incoming and outgo¬ 
ing Alfven waves, and L± is the turbulent correlation length. 
We normalize L± to a value of 75 km at the photosphere 
(Cranmer et al. 2007). This allows us to write L± = 0.15R, 
where R is the radius of the flux tube, which scales as 
Below the transition region, Q/Gphen ~ 0.2, which is simi¬ 
lar to what van Ballegooijen et al. (2011) found for closed 
loops. Above the transition region, Q/Qphen 1, which may 

be explainable if the actual correlation length L± expands less 
rapidly than we assumed from its proportionality with the flux 
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Figure 4. The time-averaged results for the PCH model are shown as a function of Alfven travel time as a proxy for height above photosphere: (a) Energy 
density, (b) rms velocity, (c) heating rate, and (d) rms magnetic field fluctuation. 


tube radius R. Because the ratio QjQ-phen strays as much as 
an order of magnitude away from unity, we further compare 
braid’s computed heating rate with analytical expressions 
that contain correcting factors that take into account efficiency 
of turbulence as a function of height. 

The first of the two efficiency factors we use is based on the 
prior work of Dmitruk et al. (2001), Dmitruk & Matthaeus 
(2003), and Cranmer et al. (2007). The extended expression, 
Qzfii is given by: 

Q'Lffl ~ t^turb^phen? (5^) 

1 

l+(feddy/hef)”’ 

where feddy is the outer-scale eddy cascade time, feddy = 
L\_ ( 1 -I-Ma) / v_L, and fref is the macroscopic Alfven wave 

reflection timescale, fref = 1/|V • Va| (Cranmer et al. 2007). 


In the expression for feddy. the velocity v_l is the ampli¬ 
tude of perpendicular fluctuations, defined previously in the 
BRAID results as Avrms- When feddy ^ het. turbulent heat¬ 
ing is quenched. The turbulent efficiency factor eturb accounts 
for regions where energy is carried away before a turbulent 
cascade can develop. The exponent n is set to 1 based on an¬ 
alytical and numerical models by Dobrowolny et al. (1980), 
Matthaeus & Zhou (1989), and Oughton et al. (2006). The 
efficiency factor works to make Qz.ov < Gphen. bringing the 
ratio QlQzfii up relative to Q/Q-phea- At low heights, where 
the efficiency factor is low, the inclusion of the efficiency fac¬ 
tor defined in Equation (5)b does help to bring Q and Qz.m 
into better agreement. At large heights, the efficiency factor 
is closer to 1, so Qphen and Qzsn both underestimate the heat¬ 
ing rate computed by BRAID. 

An alternative efficiency factor has emerged from studies 
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Figure 5. For the PCH model, we show (a) the time-averaged amplitudes of 
incoming and outgoing Alfven waves and (b) the rms amplitudes of magnetic 
and velocity fluctuations from BRAID and ZEPHYR. 


Figure 6. Similar to Figure 5; For the EQS and NAR models, we show the 
(a) outgoing and incoming Elsasser vaiiables and (b) the rms amplitudes of 
magnetic and velocity fluctuations. 


of closed coronal loops driven by slow transverse footpoint 
motions. In such models, magnetic energy is built up from the 
twisting and shearing motions of the field lines (Parker 1972), 
and the energy dissipation appears to follow a cascade-like 
sequence of quasi-steady relaxation events. Cranmer (2009) 
parametrized the time-averaged heating rate in these models 


where Ly is often dehned as the “loop length” for closed mag¬ 
netic structures, and a is an exponent that describes the sub- 
diffusive nature of the cascade in a line-tied loop. The quan¬ 
tity in parentheses is a ratio of timescales; this ratio is the 
nonlinear time over the wave travel time. To set the value of 
Lj| in our open-field models, we follow Schrijver et al. (2004), 
who found that open and closed regions can be modeled using 
a unihed empirical heating parametrization when the actual 
loop length L is replaced by an effective length scale 


L 

l + (T/To) 


(7) 


where Lq = 50 Mm.Thus, for open-field regions in which L ^ 


Lq, we use L|[ = Lq. 

The exponent a describes how interactions between 
counter-propagating Alfven wave packets can become mod¬ 
ified by MHD processes such as scale-dependent dynamic 
alignment (Boldyrev et al. 2009). The value a = Q corresponds 
to a classical hydrodynamic cascade. Gomez et al. (2000) 
constructed a model of MHD turbulence in which a = 1.5, 
and van Ballegooijen (1986) constructed a random-walk type 
model in which a = 2 (see also van Ballegooijen & Cran¬ 
mer 2008). Rappazzo et al. (2008), however, found from 
numerical simulations that a can occupy any value between 
1.5 and 2, depending on the properties of the background 
corona and wave driving. Cranmer (2009) created an ana¬ 
lytic prescription for specifying a in a way that agrees with 
the Rappazzo et al. (2008) results. Subsequently, Cranmer 
(2009) and Cranmer et al. (2013b) found that when modeling 
the coronal X-ray emission from low-mass stars, the longest 
loops—which seem to be most appropriate to compare with 
open-held regions—tend to approach the high end of the al¬ 
lowed range of exponents (i.e., a « 2). Thus, in this paper 
we use a = 2 and note that the differences between Qioop and 
Gphen will be reduced for smaller values of the exponent. At 
the largest heights, Qioop better matches the BRAID numerical 
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Figure 7. Ratio of magnetic to kinetic energy density differs from equiparti- 
tion below the transition region. 

results than either Qphen or 2 z, 07 - 

For additional comparison, in Figure 9 we plot the ratio of 
BRAID numerical results with the phenomenological heating 
rate for all three of our flux tube models. The behavior of 
the heating rate expressions with efficiency factors is similar, 
so we show only the ratio Q/Gphen in Figure 9. Note that 
the EQS model behaves similarly to the PCH model, but the 
NAR model exhibits a marked decrease in the ratio Q/Qphen in 
the low corona before increasing again to approach the EQS 
model at the top of the grid. This behavior is reminiscent of 
the closed-field models of van Ballegooijen et al. (2011), in 
which QlQ-pben came back down to values of 0.2-0.3 in the 
coronal part of the modeled loops. 

4. EXTENDED ANALYSIS OF TIME VARIABILITY 
4.1. Statistical variations as a source of multithermal plasma 

The turbulent heating simulated by BRAID was found by 
van Ballegooijen et al. (201 1) to be quite intermittent and vari¬ 
able on small scales. Eigure 10 illustrates some of this vari¬ 
ability by showing the fluctuation energy density and heating 
rate volume-averaged over the low corona (i.e., between the 
transition region at z = 0.01 Rq and an upper height of 0.5 Rq). 
This is a similar plot as Eigure 4 of van Ballegooijen et al. 
(2011). Even with this substantial degree of spatial averag¬ 
ing, the nanoflare-like burstiness generated by the turbulence 
is evident in Eigure 10. There is a large body of prior work 
concerning such intermittent aspects of turbulent heating (see, 
e.g., Einaudi et al. 1996; Hendrix & van Hoven 1996; Dmitruk 
& Gomez 1997; Dmitruk et al. 1998; Einaudi & Velli 1999). 

The time-varying heating rate should also give rise to a sim¬ 
ilarly variable coronal temperature structure. We investigate 
the possibility that the resulting stochastic distribution of tem¬ 
peratures may be partially responsible for the observational 
signatures of multithermal plasmas—e.g., nonzero widths of 
the differential emission measure (DEM) distribution. Asgari- 
Targhi & van Ballegooijen (2012) studied the spatial and tem¬ 
poral response of a conduction-dominated corona to the sim¬ 
ulated variations in Q from BRAID. They found that con¬ 



HEIGHT ([R/Rsum] - 1) 


Figure 8. We compare the time-averaged heating rate from BRAID for the 
PCH model with the phenomenological expressions given by Equations (4), 
(5) and (6) by (a) directly plotting them as a function of height and (b) plotting 
ratios of the rates, with a thin black line showing unity. 

duction leads to a “smeared out” temperature structure that 
nevertheless retains much of the bursty variability seen in the 
heating rate. Here, we perform an even simpler estimate of 
the distribution of temperatures by taking the distribution of 
volume-averaged heating rates {Q) shown in Eigure 10 and 
processing each value through the simple conductive scaling 
relation of Rosner et al. (1978). Thus, 


{f) \{Q)) 


(8) 


where {T) is an estimated volume-averaged coronal tempera¬ 
ture. The normalizing value of the heating rate (2) is assumed 
to be the mean value of {Q) seen in the BRAID simulation. 
Eor simplicity, we take the normalizing value of the tempera¬ 
ture (r) to be the maximum coronal temperature found in the 
corresponding ZEPHYR model from Cranmer et al. (2007). 
The PCH, EQS, and NAR models exhibited values of (f) of 
1.352, 1.224, and 1.675 MK, respectively. 

Eigure 1 1(a) shows the distribution of derived values of (T) 
for the PCH model. If the temperatures along this flux tube 
were measured by standard ultraviolet and X-ray diagnostics, 
with time integrations long in comparison to the scale of vari- 
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Figure 9. For each model, we show the ratio of the time-averaged heating 
rate from BRAID with the phenomenological expression given by Equations 

(4). 

ability in the BRAID model, then this distribution would be 
equivalent to the DEM. For each simulated DEM, we mea¬ 
sured its representative “width” in the same way as described 
by Schmelz et al. (2014); i.e., we used the points at which the 
DEM declined to 0.1 times its maximum value. For the PCH, 
EQS, and NAR models, we found widths of 0.2294, 0.2259, 
and 02839 in units of “dex” (log T), respectively. 

Figure 11(b) compares the properties of the three simulated 
DEMs with a selection of observationally derived coronal- 
loop DEMs from Schmelz et al. (2014). The BRAID mod¬ 
els do appear to reproduce the observed multithermal nature 
of coronal plasmas, both in the absolute values of the widths 
(which fall comfortably within the range of the observed val¬ 
ues) and in the overall trend for hotter models to have broader 
DEMs. Of course, the open-field models studied here only 
span a very limited range of central temperatures (f) in com¬ 
parison to the observed cases. 

4.2. Power spectrum of fluctuations 

For the PCH model, we investigated in detail the power 
spectrum of the velocity fluctuations caused by the Alfven 
waves. To generate the Fourier transform, we assumed con¬ 
stant time spacing using the model results spanning from 
f = ta to f = 2ta. We subtracted the mean, doubled the length 
to make a periodic sequence, and then fed the cleaned quan¬ 
tity into a traditional EFT procedure. The power spectrum is 
the product of the result of that EFT procedure with its com¬ 
plex conjugate. This procedure gives us these spectra at each 
height z. 

In Figure 12, we examine a contour plot of the power in 
these fluctuations. Certain frequencies in the 10“^ Hz to 10“' ^ 
Hz range have a relatively high amount of power at all heights, 
while at the higher frequencies, there is an increase in power 
as a function of height. This is more easily seen in Figure 
13. We show the basal input spectrum with a dashed line, 
and the power spectrum at the upper boundary of our model 
lies above the power spectrum at the photospheric base for 



0 500 1000 1500 2000 

TIME (s) 


Figure 10. Due to the random jostling that generates the Alfven waves, the 
(a) spatially-averaged energy density in the corona and (b) spatially-averaged 
heating rate per unit volume in the corona are bursty and strongly time- 
dependent in nature. 

the highest frequencies. Both boundaries show that there is a 
boost above the input spectrum for high frequencies. 

Figure 13 shows that the high-frequency part of the BRAID 
turbulent power spectrum appears to be a power law V oc. /“”, 
where the value of n varies a bit with height in the model. 
Fooking only at the FFT data with / > 0.03 Hz, we found that 
n « 4.5 at the photospheric base, and then it steepens at larger 
heights to take on values of order 5.2-6.4 (i.e., a mean value 
of 5.8 with a standard deviation of ±0.6) at chromospheric 
heights below the TR. In the corona, however, n decreases a 
bit to a mean value of 4.9 and a smaller standard deviation of 
±0.3. Some of the quoted standard deviations are likely due 
to fitting uncertainties of the inherently noisy power spectra, 
but it is clear that the corona exhibits less variation in n than 
the regions below the TR. 

Despite several decades of spacecraft observations of 
power-law frequency spectra in the turbulent solar wind, it 
was not immediately obvious that the BRAID power spectrum 
should have exhibited such power-law behavior. Spacecraft- 
frame measurements are often interpreted as being a spatial 
sample through quasi-stationary wavenumber variations (see, 
e.g., Taylor 1938; Horbury et al. 2012). In strong-held MHD 
turbulence, the dominant wavenumber cascade is expected to 
be in the k± direction. However, the Alfvenic Huctuations 
that make up an MHD cascade have a dispersion relation in 
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Figure 11. (a) Scaled histogram of the probability distribution of coronal 
temperatures estimated from BRAID heating rates. Temperatures at which 
the distribution drops to 0.1 times the maximum value (dashed line) are in¬ 
dicated by solid red lines, (b) Observed DEM widths from Schmelz et al. 
(2014), plotted as the relative width in logT space and computed from the 
reported number of logT “bins” (gray diamonds), compared with simulated 
DEM widths from the PCH (red circle), EQS (green circle), and NAR (blue 
circle) models. 



Figure 12. Contour plot of the power spectrum in the PCH model showing 
power in color, as a function of height and of frequency. 



which the frequency depends primarily on Thus, MHD 
(especially RMHD) turbulence is typically described as “low- 
frequency turbulence” and the idea of an inherent power-law 
frequency cascade is met—usually, rightly so—with skepti¬ 
cism. 

There has been one proposed model in which a power-law 
spectrum in frequency (i.e., in ky) occurs naturally and with¬ 
out the need for substantial parallel cascade: the so-called 
“critical balance” model of Goldreich & Sridhar (1995). In 
this picture, strong mixing is proposed to occur between the 
turbulent eddies (primarily moving perpendicular to the back¬ 
ground field) and Alfven wave packets (moving parallel to the 
field), such that the parameter space “filled” by a fully devel¬ 
oped cascade is determined by the critical balance parameter 


Figure 13. Power spectra for the upper and lower boundaiies of the PCH 
model compared to the input spectrum show a boost in power at high fre¬ 
quencies. 


x< 1 corresponds generally to low-frequency fluctuations 
with k|| <C as is expected in anisotropic MHD turbulence. 
This parameter can also be interpreted as a ratio of timescales. 

Following the implications of critical balance led Goldre¬ 
ich & Sridhar (1995) to a phenomenological expression of the 
time-steady inertial range, given here as 

Va v_l(A:_l) 

£(^l|,^±)= ^3 g(x) (10) 


^ ^11 Va 

k_LVx 


(9) 


taking on values x < 1 (see also Higdon 1984). The limit 


where this combines Equations (5) and (7) of Goldreich & 
Sridhar (1995). Above, £ is a three-dimensional power spec¬ 
trum that gives the magnetic energy density variance (in 
velocity-squared units) when integrated over the full volume 
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of wavenumber space, 

^ = Jd\E(^,kj^) = JdfVif), (11) 

and we also define the frequency spectrum V(f) in a similar 
way. Note that lirf = uj = Va in the local rest frame of the 
plasma under the assumption that the fluctuations are Alfven 
waves. 

In Equation (10) above, vx(A:x) is the reduced velocity 
spectrum, which specifies the magnitude of the velocity per¬ 
turbation at length scales and A:jj'. This spectrum is often 
assumed to be a power-law with vx oc k'^. The exponent 
m has been proposed to range between values of 1/3 (strong 
turbulence; Goldreich & Sridhar 1995) and 1/2 (weak turbu¬ 
lence; Galtier et al. 2000). Lastly, the function g(x) in Equa¬ 
tion (10) is a “parallel decay” function that is expected to be¬ 
come negligibly small for x 1- Because g(x) is normalized 
to unity when integrated over all x, a simple approximation 
for it is a step function. 


^(x) 


1, X<1, 
0 , X > 1 


( 12 ) 


(see also Cho et al. 2002; Cranmer & van Ballegooijen 2012). 
The above form for gix), combined with the assumption that 
the reduced spectrum vx(^_l) extends out to k± oo, leads 
to the high-frequency end of the frequency spectrum obeying 
a power law, with 

'P(f) oc /('«+l)/(m-l) ^ (13) 

The strong turbulence case {m = 1/3, 7^ oc /“^) has been stud¬ 
ied extensively both observationally and theoretically (see, 
e.g., Horbury et al. 2012). 

Of course, the BRAID models highlighted in this paper 
do not have reduced perpendicular velocity spectra that ex¬ 
tend over large ranges of k± space. The models presented 
here (similar to those of van Ballegooijen et al. (2011)) re¬ 
solve only about one order of magnitude worth of an “iner¬ 
tial range” in k± space. Eor the step-function version of g(x) 
given above, the imposition of a vx(^j_) cutoff above an ar¬ 
bitrary kmax produces a frequency spectrum Vif) that is sim¬ 
ilarly cut off above a frequency /max determined by critical 
balance (x = 1) at kx = kmax- 

In an alternative to the step-function version of gix), Cran¬ 
mer & van Ballegooijen (2003) found an analytic solution for 
six) by applying an anisotropic cascade model that obeyed 
a specific kind of advection-diffusion equation in three- 
dimensional wavenumber space. The general form of gix) is 
reminiscent of a suprathermal kappa function (see, e.g., Pier- 
rard & Lazar 2010), which is roughly Gaussian at low x and a 
power-law at large x- For x C Cranmer & van Ballegooi¬ 
jen (2003) found that gix) oc ;^'A3.s+4)/2^ where s is the ratio 
of the model’s perpendicular advection coefficient to the per¬ 
pendicular diffusion coefficient. It is still not known if MHD 
turbulence in the solar corona and solar wind exhibits a uni¬ 
versal value of s, or even whether or not s is even a physically 
meaningful parameter. Nevertheless, the wavenumber diffu¬ 
sion framework of Zhou & Matthaeus (1990) and Matthaeus 
et al. (2009) has been shown to be consistent with a value of 
s = 2 in this family of advection-diffusion equations. In a dif¬ 
ferent model of coronal turbulence, van Ballegooijen (1986) 
showed that a cascade of slow random-walk displacements of 
the held lines can be treated as the case s = 1. On the basis 


of observations alone, Cranmer & van Ballegooijen (2003) 
and Landi & Cranmer (2009) found that if s could be main¬ 
tained at small values of order 0.1-0.3, there would be suffi¬ 
cient high-frequency wave energy to heat protons and minor 
ions via ion cyclotron resonance. 

No matter the value of s or the reduced spectral index m, 
it can be shown that at large frequencies (/ > /max), a power 
law of the form gix) oc x~" produces a power-law frequency 
spectrum V oc /“" with the same exponent (see the integration 
over k± in Equation (11)). Thus, we postulate that measuring 
n from the BRAID simulations may be a way to extract infor¬ 
mation about the exponent s, with 


2 n —4 


(14) 


The values of n reported above imply a photospheric value 
of s « 1.7, which increases to s « 2.5 in the chromosphere 
(with a relatively large spread) and then decreases to s « 1.9 
in the corona. The similarities to the theoretical value of s = 2 
(from, e.g., Zhou & Matthaeus 1990; Matthaeus et al. 2009) 
are suggestive, but not conclusive. 

4.3. Nanoflare statistics of heating rate variability 

We investigated the variability of heating and energy as a 
function of height and time throughout the simulation. In a 
given finite “zone,” the energy lost via dissipative heating can 
be calculated using the heating rate as 

Eizfl) = Qiz,t)AziTTR^)At, (15) 

where the zones are defined at a set of heights, z, with un¬ 
equal spacing Az, and at a set of times, f, with equal spac¬ 
ing At = 0.25 s. In Eigure 14, we provide a contour plot of 
the energy, comparable with Eigure 6b of Asgari-Targhi et al. 
(2013). Many of the impulsive heating events that result in 
spikes of energy over a short time frame stop at the transition 
region, which lies at Q.QIRq, but some extend to ten times 
that height. There is also a lack of energy at the beginning of 
the simulation (f < ta), where the information has not yet had 
time to propagate up through the model grid. 


0.03 



1000 1500 

Simulation Time (s) 


Figure 14. Dependence of energy on height above photosphere (y-axis) and 
time (x-axis). The transition region is at 0.01 Rq. 


Eollowing the method of Asgari-Targhi et al. (2013), we 
then use box capturing to get a statistical sense of the distri- 
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bution of energy in these zones throughout the corona. In or¬ 
der to be directly comparable to their defined events, we also 
use boxes with a width of 19.4 seconds in simulation time 
and height of 19.4 seconds in Alfven travel time (a proxy 
for height). This choice in box size results initially in 118 
sections across the time dimension and 39 sections along the 
height dimension. However, the lowest 10 boxes are at heights 
below the transition region, and we plot only boxes in the 
corona in Figure 15. Additionally, we take out the first 770 
seconds corresponding to one Alfven travel time in the PCH 
model (recall Table 1 and Equation (3)) to ensure that the 
waves have had time to propagate fully throughout the corona. 

These cuts result in 78 time sections and 29 height sec¬ 
tions, giving us a total of 2262 coronal boxes. The arith¬ 
metic mean in log-space of the energy contained within these 
boxes is 21.91 with a standard deviation of 0.97. The aver¬ 
age energy contained in the boxes is lower than that found 
by Asgari-Targhi et al. (2013), since our model does not have 
two strong footpoints supplying separate sources of counter- 
propagating Alfven waves, but it is still a significant amount 
of energy. Figure 15 also shows that we find a few events 
that reach higher energies than the histogram from Asgari- 
Targhi et al. (2013), and that higher tail gets well into the 
classical nanoflare expected energies (Hannah et al. 2011). 
The minimum-energy event contains lO'®®"^ erg while the 
maximum-energy event contains lO^'*'®^ erg. It is worthy of 
note that the peak energies captured in these boxes fall well 
within the “picoflare” range (Axford et al. 1999; Parnell & 
Jupp 2000). All of these events are a natural product of the to¬ 
tal heating coming out of BRAID, suggesting that the physics 
contained within these models may lead to the formation of 
pico- and nanoflares. 

5. DISCUSSION AND CONCLUSIONS 

We have analyzed three typical open magnetic held struc¬ 
tures using one-dimensional time-steady modeling and three- 
dimensional time-dependent RMHD modeling. These struc¬ 
tures represent characteristic flux tubes anchored within a 
polar coronal hole, on the edge of an equatorial streamer, 
and neighboring a strong closed-field active region. We 
show that the time-averaged properties of the higher¬ 
dimensional BRAID models agree well with that of the less- 
computationally-expensive ZEPHYR models. 

We looked in detail at the energy partitioning, as BRAID 
imposes no assumptions or restrictions. At heights above 
the transition region, equipartition is shown to work well 
to describe the results from BRAID, and is assumed in the 
ZEPHYR algorithm. We also compared the energy partition¬ 
ing with predictions from non-WKB reflection for a range 
of frequencies, showing the limitations of the linear method 
for such predictions. The time-averaged heating rates from 
BRAID are lower than the phenomenological expression in 
Equation (5) for the heating rate below the transition region, 
rising sharply toward the upper boundary. 

In BRAID, Alfven waves are generated by random foot- 
point motions, whose properties and driving modes are de¬ 
scribed in Section 2.2. As the Alfven waves propagate up¬ 
ward from the photosphere to the open upper boundary at a 
height of 2 solar radii, they partially reflect and cause tur¬ 
bulence to develop. With the time-dependence included in 
BRAID, we are able to show the bursty nature of turbulent 
heating by Alfven waves. We show that this heating brings 
energy up into the corona and provide a statistical distribution 
of energy per event. The more energetic events (i.e., boxes 


with the most energy) fall within expected nanoflare values. 

Overall, we show that time-steady modeling does a good 
job of predicting the time-averaged results from time- 
dependent modeling. There is, however, a bounty of infor¬ 
mation that can be found only by looking at changes in the 
heating rate over time. Moving from one dimension to three 
allows the model to contain more realistic physics. We have 
shown that these models of typical magnetic held structures 
provide additional compelling evidence to support the idea 
that Alfven-wave-driven turbulence heats the corona and ac¬ 
celerates the solar wind. 
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